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Abstract. A one-dimensional rule-based model for flocking, that combines velocity 
alignment and long-range centering interactions, is presented and studied. The induced 
cohesion in the collective motion of the self-propelled agents leads to a unique group 
behaviour that contrasts with previous studies. Our results show that the largest 
cluster of particles, in the condensed states, develops a mean velocity slower than the 
preferred one in the absence of noise. For strong noise, the system also develops a 
non-vanishing mean velocity, alternating its direction of motion stochastically. This 
allows us to address the directional switching phenomenon. The effects of different 
sources of stochasticity on the system are also discussed. 

PACS numbers: 05.40.-a, 05.70.Fh, 64.60.Cn, 87.10.-e 
1. Introduction 

The phenomenon of flocking (also known as swarming) refers to the emergent collective 
motion regularly seen in groups of animals such as bird flocks, fish schools, bison or sheep 
herds, insect or bacteria swarms, etc. [Il[2l|3llll[5l[6l[7] In simple terms, the formation 
of these sort of groups in nature, implies some kind of condensation in position and/or 
velocity spaces. The mechanisms or behaviours followed by the individuals that compose 
the group, that lead to this type of condensation, have been classified by Reynolds [S] 
as: the desire to match the velocity of flockmates (alignment), the desire to stay close to 
flockmates (centering), and the desire to avoid collisions (separation). These definitions 
have been widely adopted by researchers from a variety of disciplines. 

In the realm of the physical sciences, the study of flocking has greatly benefited 
from the insight of statistical mechanics to describe dynamical phase transitions and 
cooperative phenomena p]. The interest in the study of this kind of systems comes from 
the fact that they represent prototypical out-of-equilibrium systems. The formulations 
to study the flocking phenomenon can be succinctly classified in rule-based [TOl [HI [121 
tta [H [la [161 El HE], Lagrangian (trajectory-based) |l9l|20l|2ll[22l[23l[2a[25l|26] 
and Eulerian (continuum) models f27\ [2S1 EH]- Regarding the dimensionality of the 
models developed, most of them have been defined in dimensions higher than one 
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[IQllIIlliaiTSlITlIISllMlEIlESlEalEaE^ since the velocity of the 

self-propelled particles (SPP) in these models can have continuous values. In contrast, 
only a few one-dimensional (ID) models have been studied [5l [161 113 UHl [191 [26l|3T]. For 
these, the direction of the particles velocity can only take discrete values, in fact, only 
two. Nonetheless, recent experimental and theoretical studies by Buhl et al [S] and Yates 
et al [6| have proven the usefulness of ID models in elucidating the directional switching 
phenomenon. This phenomenon regards a very interesting aspect of the movement of 
many animal groups, when they suddenly switch their direction of motion in a coherent 
form. This seems to be an intrinsic property of the collective motion of animals as it 
has been observed, e.g., in groups of locusts and starlings (see |H] for more details). 

Among the one-dimensional models, one can find of the Lagrangian kind like 
the one by Mikhailov and Zanette [12], that describes a population of self-propelled 
particles (SPP) with attractive long-range interactions, i.e., includes centering as the 
only ingredient from the Reynolds considerations. Here, the particles try to attain 
one of two preferred velocities [vq or -vq) through a non-linear friction term cubic in 
nature. The system presents multistable dynamics and can be found either in coherent 
traveling motion {ordered state) where all the particles move together in a condensed 
group, developing a non-zero mean velocity with the same direction, or in an incoherent 
oscillatory state {disordered state) where the particles perform random-phase oscillations 
around a fixed position with zero mean velocity. By increasing the additive noise 
intensity, the system can be driven from the coherent traveling state to the oscillatory 
state through a discontinuous phase transition. However, for noise intensities below the 
critical point, the two states are accessible depending on the initial conditions since the 
system shows a strong hysteresis. 

On the other hand, one can find rule-based models like the ID version of the 
famous Vicsek et al model [TU], by Czirok et al [21], that describes a population of 
SPP that interact only through condensation in velocity space (alignment according to 
Reynolds). Via short-range interactions and a majority rule, the particles try to match 
the velocity of their nearest fiockmates. The system exhibits spontaneous symmetry 
breaking and self-organization depending on the intensity of the also additive noise, 
leading to a continuous kinetic phase transition from a homogenous to a condensed 
phase. In contrast with the 2D Vicsek model [ID] where the speed of the particles 
is fixed, i.e., |fo| = constant, in the ID version [31], at each time step, particles try to 
acquire a velocity close to the two prescribed ones {vq or ~vq) through an antisymmetric 
function G{u) that depends linearly on the local mean velocity {u{t))^ of the neighbours 
of a given particle i. This function can also be written in a continuous form with no 
change in the scaling properties of the system as claimed in [31] . 

Let us also consider the lattice model by O'Loan and Evans [T^] where the 
velocity of the particles can take three discrete values (—1,0,1). Here, through an 
alignment interaction in the form of a majority rule, that is applied asynchronously 
with probabilistic updates (particles are selected and updated randomly in a Monte 
Carlo fashion), the system shows a continuous phase transition from a homogeneous to 
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a condensed phase similar to that of the previous modeL However, and in contrast with 
the previous model, the condensed phase is not symmetry-broken but alternating. We 
must advert the two sources of noise in this model — coming from the way particles 
are updated and how the alignment rule is applied — since no additive noise is present. 
The later generalization of this model by Raymond and Evans [T7| includes all three 
Reynolds' flocking behaviours that translates in richer dynamics: the system also shows 
an homogeneous flock with homogeneous density and a fixed global non-zero velocity, 
and dipole states. However, the different interactions present in these models are local, 
i.e., short in range regarding the spatial coordinate. 

Finally, it is worth to point out the variations of the two previous rule-based 
models introduced by Buhl et al [5], Yates et al |6] and Bode et al [18], that also show 
an alternating condensed phase that resembles the directional switching phenomenon 
mentioned before. In the first. Buhl et al consider a generalization of the Czirok et al 
model, that weights the local mean direction of motion of the neighbours of a given 
particle with that of the particle itself in a deterministic way. In the second, Yates 
et al adapted the Buhl et al model with an additive noise whose intensity is weighted 
with a non-trivial function of the local mean velocity: the individuals increase the 
randomness of their motion when the local alignment is low. In the third. Bode et 
al mix features from the Czirok et al , O'Loan and Evans, and Raymond and Evans 
models, concentrating on how noise enters in the system. There, they have identified 
three different approaches (that can be generalized to higher dimensions) to include 
stochastic errors in ID SSP models and classified them as: 

(z) Adding a random variable to the preferred direction of individuals (i.e., the 
inclusion of additive noise as in [THIEI]). 

{ii) Asynchronous and probabilistic updates (as in [161 113 )• 

{Hi) Varying the probability and accuracy with which individuals execute their 
behavioural rules (as in [THl [H] also). 

In this way. Bode et al have included local alignment and centering interactions in 
their model, that are applied taking inspiration from the approaches [ii) and {Hi). 
Nonetheless, the interactions considered in all of these models are short-ranged from 
the point of view of the position of the particles as well. 

Our aim in this work is to extend the ID model of Czirok et al [21] by providing 
it with a long-range centering interaction (much in the spirit of Mikhailov and Zanette) 
along with the already included local alignment interaction and additive noise. We 
demonstrate a variety of new flocking regimes: some of them resemble those found in 
the models described in the previous paragraph. We also show that the noise induced 
phase transition, from the homogeneous to the symmetry-broken condensed phases, 
changes into what seems a phase transition between two kinds of condensed states: one 
with broken symmetry and one with an alternating flock that moves with non-zero mean 
velocity, changing its direction of motion stochastically. The latter seems stable even for 
strong noise and different system sizes. Additionally, we explore two ways to implement 
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Figure 1. In the lower plots, the dynamics of the ID Czirok et al model for L = 1000, 
N — 1000 and pa = 1, for mcreasing values of 77 (from left to right 77 = 0,3,5,7, 
respectively) . The darker grey level represents higher particle density and the systems 
were let to evolve 15000 time steps. The upper plots present the corresponding 
normalized distribution P{u) of the velocities Ui of the individual particles in the 
steady state. Notice how the steady state of the system changes from one with broken 
symmetry to a homogeneous one as the noise amplitude increases. 



the centering and alignment interactions — one probabilistic and one deterministic — 
that allows us to address the relation of the directional switching phenomenon with 
different sources of stochasticity in ID SPP models. 

2. Model 

The ID model by Czirok et al consists of N off-latice SPP along a line of length 
L. The particles are characterized by their coordinate Xi and dimensionless velocity Ui 
updated as 

Xi{t + 1) = Xi{t) +VoUi{t), 

u,{t + l) = G{{u{t))^)+^i. 

The local average velocity {u{t))- for the ith particle is calculated within an interval 
[xi — A,Xi + A], where A = 1 is fixed, and updated for all particles at every time step. 
Propulsion and friction forces are incorporated through the antisymmetric function 
G which sets the velocity in average to a prescribed value Vq, in such a way that 
G{u) > M for < M < 1 and G{u) < m for m > 1. In this case, one cannot apply 
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the constant velocity constraint of [lOJ since it strongly discretizes the system, leading 
to a discrete state cellular automaton. It is also important to remark that, similarly with 
the Raymond and Evans model [TTj, function G{u) also works as the interaction term. 
The noise, is distributed uniformly in the interval [—ri/2, r]/2] and affects directly the 
velocity of the particles, while an update time At = 1 has already been applied. 

In this work, t>o is kept constant {vq = 0.1), and the adjustable control parameters 
for the model described above are the average density of the particles, pa = N/L, and 
the noise amplitude t]. Function G was implemented in [31] in a very simple way as, 

G{u) = l[u + sign(M)] . (2) 

The results do not depend on the properties of G{u) for m as stated in [3l]. In our 
case, simulations with continuous choices of G{u) were carried out leading to similar 
results. Random initial conditions are applied along with periodic boundary conditions. 

In the lower plots of figure [T| we show the time evolution of the density of particles, 
p{x,t), and the normalized distribution P{u) of the velocities Ui of the particles, in the 
upper plots, for the model given in ([T| and ^ with L = 1000, = 1000 and increasing 
values of the noise amplitude 77(6 [0, 7]). The corresponding average density of particles 
is Pa = 1- Depending on the chosen noise amplitude, 77, the system can be driven through 
an order- disorder continuous phase transition. For example, for rj < rjc ^ 5, the system 
reaches an ordered state in a short time (see figures [ij^a) and[l](b)), characterized by 
a spontaneous symmetry breaking and clustering of the particles. On the contrary, for 
larger values of rj {> rjc), a disordered phase can be found characterized by a random 
velocity field: see, e.g., figure [l](d). As one increases r/, one can appreciate how the 
distribution P{u) in the steady state, that starts single-peaked and narrow for 77 = 
(see upper plots in figure [T]), flattens and broadens out when going from the ordered to 
the disordered state (from left to right). 

In order to quantify the time evolution of the system, we define the instantaneous 
order parameter as 



N 



j:u,it) 



=1 



(3) 



and 



We also measured the instantaneous position and velocity of the centre of mass, given 
by 

1 ^ 

xcuit) = -Y^^iit)^ (4) 
1=1 

1 ^ 

ucuit) = —J2ui{t), (5) 
1=1 

respectively. 

Before we proceed, there is one thing that needs to be clarified. Due to periodic 
boundary conditions, a ID system has the topology of a ring. This implies there are 
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Figure 2. In (a), (b), (c), (e) and (f) the steady-state accumulated order parameter 
(f>ss: given in (|6]), is plotted as a function of the noise amplitude rj. In each case, the 
system was let to evolve over 50000 time steps until it reached the steady state, and 
the order parameter was accumulated over 250000 time steps after the transient. In (a) 
the ID Czirok et al model, given in ([T]) in its original form, is considered for different 
values of Pa and N — 1000. In (b) and (c), and the details (e) and (f) corresponding 
to the grayed out portions of (b) and (c), respectively. Version 1 (VI) and Version 2 
(V2) of our model, given in ([7| and ([s]), are considered for = 1, = 0.05,0.1 and 
different system sizes {N — 500, 1000,4000). In (d) the instantaneous velocity of the 
centre of mass ucuit), given in ([5|, is plotted for pa = 1, /i = 0.1, rj ~ 7 and different 
system sizes in the steady state: the black solid lines correspond to Version 1 of our 
model while the grey solid lines to Version 2 (see the text for more details). 



two distances (running in opposite directions) from any particle to any specified point 
in the ring. Additionally, let us not forget that random initial conditions are taken both 
in the position and in the initial direction of the particles velocity, and that alignment 
interactions favour the formation of clusters as shown in figure [Tj For simplicity, the 
position of the centre of mass, as given in (|4]), is determined with the particles position, 
Xi{t), measured with respect to the origin from which the length L of the system is 
measured as well. However, this does not affect the overall results obtained as we will 
see below. 

To account quantitatively for the transition from the ordered to the disordered 
state, in figure we plot the steady-state accumulated order parameter, defined as 

1 f'^ 

'^^^ = T-SL f Jo ^^^"^ 
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Figure 3. Schematic diagram of the small velocity regime compared to the large 
velocity regime. 

as a function of the noise amplitude rj for = 0.5, 1, 2. This figure makes evident the 
dependence of the critical point rjc on the average density of the system (see more 
details in [HT]). 

Let us now provide the ID Czirok et al model with a new interaction rule that 
considers centering along with alignment. However, there are some constraints that 
have to be taken into account. The dynamics in the Czirok et al model may be sensitive 
to the prescribed velocity fo of the particles. As reported in [12] for the original two- 
dimensional Vicsek et al model [TU] — also known as the scalar noise model (SNM) 
— an important quantity to be considered is the ratio of the interaction radius and 
the velocity of the particle, R/vq, compared to the update time At in the numerical 
simulations. The SNM model was proposed to study the motion of bird flocks and/or 
bacterial colonies. The motion of the particles in such systems is quasi-continuous, i.e., 
the reaction time of the birds is, usually, significantly faster than the characteristic time 
that is needed to travel through their interaction radius R. This condition imposes the 
following constraint: At ^ R/vq- In our case, once fixed the interaction radius A = 1 
and the update time At = 1, the above condition becomes fo <^ 1, the same as in the 
SNM model. This velocity domain is referred as the small velocity regime (see figure 
[3|. If the velocity of the particles is allowed to become larger (t>o > 0.3, referred as 
the large velocity regime in [12]), the SNM model shows density waves that appear due 
to an emergent influence of the periodic boundary conditions. Moreover, for very large 
particle velocites (e.g., fo = 10) the order-disorder transition exhibits a discontinuous 
order parameter characteristic of a first order phase transition (see [12] )• 

As far as we know, the properties of the Czirok et al model, as given in Q and 
(|2|, have not yet been studied in detail outside the small velocity regime. Nonetheless, 
it is known that for certain combination of the parameters, this model can develop 
directional switching (see, e.g., [T8]). 

Taking into account the previous statements, we will keep the system in the small 
velocity regime {vq <^ 1), as we add centering to the type of interactions experienced 
by the particles. In this way, on top of the majority rule that allows particles to match 
the velocity of their neighbours (interacting in the velocity space), we will introduce a 
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simple attractive long-range rule (an interaction in the position space) that resembles 
that one used in the Mikhailov-Zanette model pi9]. However, in contrast with the 
pairwise linear forces that attach the particles in the Mikhailov-Zanette model, in order 
to determine the new direction for a given particle, our centering-rule for the Czirok 
et al model will only take into account the shortest distance from the position of that 
particle relative to the position of the centre of mass of the group xcu] this is done so 
because of the ring topology of the system. Additionally, to make the implementation 
of this centering-rule simple, the mean velocity of the particle's neighbours (obtained 
through function G{{u{t))-)) will then be directed towards the centre of mass whenever 
this rule is applied: we will call this function H{{u{t))- ;xcm)- In this way, we will be 
able to keep the velocity of the particles within the small velocity regime, so they will 
not be allowed to abandon their interaction regions in less than the characteristic time 
defined by the quantities A, fo, and At. This will also allow us to define two versions 
of the Czirok et al model, that include centering, in which these two rules are applied 
taking inspiration from two of the approaches to include stochastic errors as classified 
by Bode et al [18], listed in the Introduction of this work. 

In Version 1 (VI) of our model we consider a probabilistic approach to apply 
the alignment and centering rules — similar to approach (ii) according to Bode et 
al (see above) — but with synchronous updates. To weight one rule above the other we 
introduce the centering parameter /i in the following way: 

• With probability (1 — /i), a given bird in the flock will follow the original majority 
rule G{{u(t))-), trying to match its velocity with the mean velocity of its neighbours 
(itself included). 

• With probability fi, it will acquire the magnitude of the mean velocity of its 
neighbours, but will move towards the centre of mass of the flock instead, following 
rule H{{u{t))- ; xcm) as described above. 

We also keep the additive noise as in the original Czirok et al model. Thus, in this 
version of our model, the coordinate and dimensionless velocity Ui of the particles are 
updated as 



where / is replaced by rule G{{u(t))-) or H{{u(t)). ; xcm) according to the value of /i in 
an probabilistic way. 

In Version 2 (V2) of our model only the additive noise C,i is kept as a source of 
stochasticity, applying the rules in a deterministic fashion like in the model of Buhl et 
a/ [5] . In consequence, in this version of our model, the coordinate Xi and dimensionless 
velocity Ui of the particles are updated as 



Xi(t + 1) = Xi(t) + VoUiit), 
Uiit + 1)= 



(7) 



Xi(t + 1) = Xi(t) + VQUi(t), 

Ui{t + !) = (!- fi)G{{u{t)),) + fiH{{u{t)), ; Xcm) + ^i- 



(8) 
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Figure 4. The dynamics of Version 1 of our model, given in ([7|, is presented in the 
lower plots for L = 1000, N = 1000, pa = 1 and /i = 0.05, for increasing values 
of rj (from left to right rj = 0,3,5,7, respectively). The darker grey level represents 
higher particle density and the systems were let to evolve 15000 time steps. The 
upper plots present the corresponding normalized distribution P{u) of the velocities 
Ui of the individual particles in the steady state. Please notice how, in the absence of 
noise (77 = 0), P{u) shows two peaks centered (from left to right) around the velocity 
of particles not yet absorbed by the densest cluster and around the velocity of the 
densest cluster, respectively (see the text for more details). 



In this way, in both versions of our model one can go from the original Czirok et al model 
(for /i = 0) to one that resembles the interaction type present in the Mikhailov-Zanette 
(for /i > 0). This allows to continuously go from a purely alignment type of model to 
a purely centering one (for /i = 1). However, let us not forget that the scaling in the 
Czirok et al model depends on the average density of the particles pa and the size of 
the system L, characteristic not shared by the Mikhailov-Zanette model where none of 
its properties depend on the density of particles. 

3. Results and discussion 

We will now consider the effects of centering and alignment in the Czirok et al model by 



studying cases with p > 0, first for Version 1 of our model in subsection 3.1, Later we 



show in subsection |3. 2| that , aside from some minor differences. Version 1 and Version 2 of 
our model are qualitatively equivalent regarding the collective behaviour of the system. 



We also contrast our results with those for some other models. Finally, in subsection 3.3 
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we present the scaling properties for both versions of our modeL 
3.1. Condensed states and directional switching 

In the lower plots of figure |4| we show the time evolution (i.e., the evolution of p{x, t)) of 
Version 1 of our model given in ([T]), for L = 1000, = 1000 and /i = 0.05, and increasing 
values of the noise amplitude ?7(g [0, 7]) from right to left, respectively. The symmetry 
breaking can still be appreciated as shown in the lower plots of figures |4]( a) and|4|^b). 
By increasing the noise amplitude, rj, the system is driven from a symmetry-broken 
ordered phase to an alternating flock phase, shown in the lower plots of figures ^c) 
and|4]^d). It is very clear the effect of the centering rule, eventually gathering all of the 
particles in one big cluster. Moreover, one can appreciate how the cluster with highest 
density (darkest in the figures), and thus always close to the centre of mass, develops a 
slower mean velocity with respect to the rest of the clusters: a close inspection of the 
slope of this cluster accounts for it, with its slope being higher compared to the slope of 
other clusters. This can be better appreciated by looking at the distribution P{u) of the 
velocities of the individual particles, shown in the upper plots of figure |4j In particular, 
notice the bimodal distribution of figure |4](a), where the narrow peak centered at — 1 
corresponds to those particles not yet absorbed by the densest cluster, while the wider 
one, centered a little to the left of —0.9, corresponds to the latter. As a result, the order 
parameter in the steady state, (pss, never reaches a value equal to 1, even in the absence 
of noise [rj = 0). On the other hand, as the noise intensity increases, 0ss never vanishes 
as the system develops a mean velocity with a positive magnitude in the alternating 
flock phase and contrary to what happens in the absence of centering. See, e.g., figure 
for Pa = 1 and /i = 0.05 where different system sizes are considered (lines with 
solid symbols). 

These effects are only enhanced as fi increases. For n = 0.1 and rj = 0, the order 
parameter, 0ss, reaches a smaller value in comparison with the case n = 0.05. In contrast, 
for strong values of the noise intensity rj, the alternating flock phase develops an order 
parameter with higher values than for p = 0.05. See the lines with solid symbols in 
figure [2]^c) (pa "was kept constant in both cases). 

To better illustrate this, the lower plots in figure |5] show the dynamics of Version 
1 of our model with p = 0.1. As expected, as p increases, both condensed phases are 
reached faster. Yet again, the densest cluster develops a slower mean velocity than the 
"free" particles that have not yet been absorbed by it, shown by the double-peaked 
distribution in the upper figure [Sj^a), with a wide peak centered a little less than 0.9 
and a narrow one centered at 1, respectively. As noise increases, for strong values of 
the noise intensity rj, the alternating fiock phase becomes more evident as shown in the 
lower plots of figures |5](c) and |5](d). As apparent from the figures, the frequency of 
the changes in direction of the densest-cluster mean velocity increases with rj, while the 
distribution P{u) broadens out (see the upper plots in figure [s]). 

In order to explain this, let us consider the combined effects of the alignment and 
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Figure 5. The dynamics of Version 1 of our model, given in ([7]), is presented in 
the lower plots for L = 1000, N = 1000, Pa = 1 and /i = 0.1, for increasing values 
of rj (from left to right rj = 0,3,5,7, respectively). The darker grey level represents 
higher particle density and the systems were let to evolve 15000 time steps. The upper 
plots present the corresponding normalized distribution P{u) of the velocities Ui of the 
individual particles in the steady state. As rj increases, the frequency of changes in the 
direction of motion of the densest cluster increases as well. 

centering interactions on the densest cluster, first for tlie case witliout noise [r] = 0). 
Because of tfie presence of centering, particles in this cluster consistently move towards 
the centre of mass even against the direction of motion of the majority. This results in 
a slowdown of the mean velocity of the group j|] See, e.g., the steady state of the case 
fi = 0.1, where an ordered phase with broken symmetry develops, moving to the left as 
shown in figure |6](a). For particles at the front of the densest cluster, the centre of mass 
lies behind as illustrated with the white solid line in the detail given in figure [6|^b). In 
consequence, particles at the front turn back to the right from time to time, opposing 
the direction of the majority. However, the alignment interaction keeps particles moving 
to the left most of the time, switching the direction of the opposing particles once again. 
This means that any given particle moves mainly in the direction of motion of the group, 
but it also moves in the opposite one as corroborated by looking at P{u) in figure [6|^c), 
depicted with the solid black curve for this case. There, the higher peak is centered 
around the emergent mean velocity of the group, while the shorter one (in a symmetric 
position) appears due to the opposing particles. 

I Please remember that the alignment interaction determines the instantaneous velocity of a given 
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Figure 6. (a) The dynamics of Version 1 of our model, given in ([7]), is presented for 
L = 1000, N — 1000, pn — 1, fi — 0.1 and ij — 0. The darker grey level represents 
higher particle density. The system was let to evolve 25000 time steps, while the 
steady state measurements were taken above the horizontal dotted line corresponding 
to t = 20000. (b) Detail of the dynamics, marked with a dashed square in (a). The 
trajectory of the centre of mass, xqm, is depicted with a white solid line, (c) Semi-log 
plot of the normalized distribution P{u) of the velocities of the individual particles, 
in the steady state, for Version 1 (solid black) and Version 2 (solid grey) of our model 
with the same parameters as in (a); Version 2 of our model is given in (11 . The 
peaks of P{u) for each version of our model are centered differently (symmetric vs. 
non-symmetric) due to the way alignment and centering interactions are implemented 
(see the text for more details), (d) Instantaneous order parameter (j){t), given in (|3|, 
for Version 1 (solid black curve) and Version 2 (solid grey curve) of our model. The 
horizontal white and black dashed lines correspond to the order parameter 4>ss, for 
Version 1 and Version 2 of our model, respectively, accumulated in the steady state 
with the same parameters as in (a). Due to the extra source of stochasticity, in the 
probabilistic way alignment and centering interactions are implemented in Version 1 
of our model, the fluctuations of the instantaneous order parameter (f>(t) are stronger 
than in the deterministic Version 2. 



In this way, once a cluster is formed around the centre of mass, this will grow denser 
over the rest of the systems' time evolution, gathering more particles together and, at 
the same, hindering their abihty to fully develop the preferred velocity as dictated by 
function G{u), given in ([2]). In some sense, inside the densest cluster, the combination 
of centering and alignment works as a source of "noise" itself (also established in [T8]). 
inducing fluctuations in the steady-state mean velocity of the group even for rj = 0. See, 
e.g., figure |6](d) where the instantaneous order parameter 0(t) is plotted as a function 
of time (solid black curve). The horizontal white dashed line corresponds to the value 

particle from the local mean velocity of its neighbours through function G{u) as G{{u{t))^). 
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Figure 7. (a) The dynamics of Version 1 of our model, given in ([t]), is presented for 
L = 1000, N — 1000, pn — 1, fi — 0.1 and ij — 7. The darker grey level represents 
higher particle density. The system was let to evolve 25000 time steps, while the 
steady state measurements were taken above the horizontal dotted line corresponding 
to t = 15000. (b) Detail of the dynamics, marked with a dashed square in (a). The 
trajectory of the centre of mass, xcm, is depicted with a white solid line. Notice how, 
on one side, the densest cluster effectively oscillates around the center of mass, while 
a tail of defecting particles counterbalances its mass on the other, (c) Semi-log plot 
of the normalized distribution P{u) of the velocities of the individual particles in the 
steady state, (d) Instantaneous order parameter given in ([s]), for the whole time 
evolution of the system, (e) Instantaneous velocity of the centre of mass ucmit), given 
in ([5|, for the steady state of the system. In (d) and (e), the horizontal grey lines 
corresponds to the accumulated mean value (/)ss, and ±0ss, respectively. 



of (f)ss measured in the steady state, above the horizontal dashed hne shown in figure 
|6]^a). Regarding the "free" particles not yet absorbed by the densest cluster, in the early 
stages of evolution, the effect of centering on their speed is negligible compared to the 
alignment interaction, due to the much lower density in these regions. In the end, the 
magnitude of the emergent velocity of the densest cluster, when the noise is weak or 
absent, will depend on the balance of the centering and alignment interactions. 

Let us now consider the case when the noise is strong. Here again, the combined 
effects of alignment and centering are stronger on the densest cluster. In figure [7](a), 
the evolution of p{x,t) is shown for a system with fi = 0.1. One can appreciate that, as 
soon as a large and dense cluster is formed, it starts to change its direction of motion 
with some periodicity. In consequence, the system develops an alternating fiock phase 
that resembles that in the models of Evans et al [161 [13 • To quantify this, it is better 
to look at the instantaneous velocity of the centre of mass, mcm, given in ([s]). As shown 
in figure l7](e), in the steady state, the velocity of the centre of mass develops a constant 
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Figure 8. Semi-log plots of the density of particles p, as a function of the position x, 
taken over 62 consecutive times steps. In (a), (b) and (c), three cases are plotted 
before, at, and after the turn marked with a small arrow in figure [Tf^b), around 
t = 18475, 18847, 19095, respectively. The levels of grey correspond to the local mean 
velocity, uioc, as shown in the scale at the bottom- right, while the horizontal black and 
white arrows point in its direction. Notice the development of domains with different 
local mean velocity. 



mean magnitude (see the horizontal grey hnes) but with an alternating direction. In 
contrast with the zero-noise case, due to the presence of noise [rj > 0), particles are able 
to drift away from the centre of mass in a "diffusive" fashion. This can be observed 
when looking at the steady state of the system, above the horizontal dotted line in 
figure [7](a). As more particles leave the densest cluster, this becomes dilute and wider 
as the centre of mass is "pulled" to its back as shown by the solid white line in the 
detail plotted in figure ^h). The turn marked with the small arrow in this figure is 
analyzed in more detail in figure [s] There, the density of particles (accumulated over 
62 consecutive time steps) is plotted as a function of the position around three different 
configurations: before, at, and after the turn. It is clear the development of domains 
with different local mean velocity, uioc, corresponding to the different levels of grey. 
As shown in figures [sj^a) and[8|^c), before and after the turn, when the centre of mass 
is closer to the densest cluster (see figure [7](b)), the latter is surrounded by domains 
with small mioc- However, at the turn, the densest cluster is at its farthest position 
relative to the position of the centre of mass (see figure [7|^b)). A fluctuation develops 
at the front of the flock, followed by a small-wioc domain as shown in figure |8]^b). This 
fluctuation moves against the direction of the flock, gaining density and momentum, 
finally inverting the direction the flock as it passes through. This mechanism effectively 
renders the particles in the system oscillating around the centre of mass, with the densest 
cluster on one side, and the defecting particles counterbalancing its mass on the other 
as can be appreciated in figure ^h). In this way, the alternating flock phase in our 
model resembles the oscillatory phase in the Mikhailov and Zanette model [19], that 
also considers long-range centering but not alignment. 

Two conclusions can be drawn from this analysis. The first relates the frequency 
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of the turns with the intensity of the noise. As explained before, the presence of noise 
allows particles to drift away from the densest cluster that always remains close to 
the centre of mass. The more intense the noise is, the faster particles are able to move 
away. Consequently, the frequency of the turns should increase with the noise amphtude 
as figures [sj^c) and [sj^d) can confirm, a trend that is consistent with previous results; 
see, e.g., [HI [161 IE]- The second regards the development of a mean velocity different 
from zero by the densest cluster when strong noise is present. In contrast with the 
case without centering (/i = 0), where strong noise makes the mean velocity of the 
flock vanish in a continuous phase transition, when centering is introduced (/i > 0) the 
accumulated order parameter is far from zero - it even increases with /i! 

For strong noise, centering in combination with the alignment interaction provokes 
an ordering effect on the densest cluster. This can be understood considering that, due 
to centering, particles try to move towards the centre of mass of the group eventually 
gathering together. With the aid of the alignment interaction, correlations in the 
velocities Ui of the particles develop that cannot be destroyed by simple additive noise. 
This is due to the long range character of the centering interaction and the high density 
of this cluster. In contrast, the alignment interaction being short-ranged in nature, is 
more sensitive to the noise when taken alone. This may also explain the increase of 0ss 
with /i, being centering the dominant interaction for strong noise. 



3.2. Different sources of stochasticity and free flocks 

Regarding the collective behaviour of the system. Version 2 of our model shows the 
same overall features as Version 1. For that matter, the dynamics of four systems with 
different combinations of the parameters A^, pa and /i, is presented in figure |9] for both 
versions, side by side, in two limit cases: rj = and rj = 7. These examples illustrate 
the qualitative equivalence (and, in great measure, quantitative too) between the two 
versions of our model — as apparent from the figures — including the alternating flock 
phase. The case = 1000, Pa = 1 and p = 0.4 in the third row, is of particular interest 
since this phase is absent even for strong noise {rj = 7). We will discuss this case later. 



More details of the dynamics of Version 2 are presented in figure 10 for /i = 0.15, 
with the same values of pa, A^, L and rj as in the two cases presented for Version 1 in 
figures |4] and |5} Again, the system shows collective behavior equivalent to that of Version 
1. Moreover, the mechanism that produces the alternating flock phase in Version 2, is 



exactly the same as the one described in subsection Nonetheless, there are some 
differences that need to be pointed out. The first regards the distribution P{u) of the 
velocities of the individual particles for the case rj = 0, shown in the upper plot of figure 



10 a), as this exhibits two peaks centered around —0.9 and — 0.6|§] As in Version 1, in 
Version 2 of our model, particles moving in the direction of the densest cluster move 
with speeds close to |1|, in contrast, particles moving in the opposite direction move 



§ A third peak should appear around —1, corresponding to particles not-yet-absorbed by the densest 
cluster. This peak is not present here since there were no such particles when P{u) was measured. 
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Figure 9. Dynamics of Version 1 and Version 2 of our model for different combinations 
of the parameters N, and /x in two limit cases: 77 = and 77 = 7. The darker grey 
level represents higher particle density. The case on the third row is of particular 
interest since the alternating flock phase is fully suppressed even for strong noise (see 
the text for more details). 



with speeds around |1| — 2//. These values correspond to the centers of the two peaks 
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Figure 10. The dynamics of Version 2 of our model, given in ([8]), is presented in 
the lower plots for L = 1000, N — 1000, Pa = 1 and fi = 0.15, for increasing values 
of rj (from left to right rj = 0,3,5,7, respectively). The darker grey level represents 
higher particle density and the systems were let to evolve 15000 time steps. The upper 
plots present the corresponding normalized distribution P{u) of the velocities Ui of 
the individual particles in the steady state. Aside from the form of the distribution of 
P{u), this version of our model shows the same overall features as Version 1. 



of P{u) shown in figure [lO|(a). This comes as a resuh of the deterministic way the 
ahgnment and centering rules are apphed in Version 2; see ([s]). Also, please notice that 
these peaks are narrower and comparable in height, in contrast to the corresponding 
ones for Version 1 — wider, symmetric in position, but highly asymmetric in height — 
as illustrated in figure ^c) for both versions with the solid grey and black filled curves, 
respectively (/x = 0.1 in this case). 

This comes about the fact that, in both versions of our model, the combination 
of alignment and centering is a source of noise in itself — all the more so, considering 
the additive noise is absent for rj = 0. However, Version 1 of our model has an extra 
source of stochasticity in the probabilistic way these interactions are applied; see ([T]). 
As a result, the fluctuations of the instantaneous order parameter 0(t) are stronger in 
Version 1 than in Version 2, as shown in figure [6]^d) for /i = 0.1 with the solid black and 
grey curves, respectively. Before we go on with the discussion, let us provide a context 
regarding the development of directional switching and its relation with different sources 
of stochasticity in ID SPP models. 

In [18j, Bode et al state that the addition of simple noise terms in the particles 
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direction of motion is not necessarily sufficient to describe and explain collective motion 
in animal groups. This remark is made in regard of the directional switching shown by 
the marching locusts in the experiments by Yates et al reported in |6j. In the latter, a 
variation of the Czirok et al model is proposed to explain this behaviour, where the local 
alignment interaction competes with each particle's own velocity in addition to a non- 
trivial additive noise. In response. Bode et al proposed a fully stochastic ID SPP model 
that does not consider any kind of additive noise, but an intrinsic source of stochasticity 
coming from an asynchronous updating with random sampling of neighbours, along 
with local concentric centering and alignment interaction zones. From both versions of 
our model, one can appreciate that the stochasticity coming from the combination of 
alignment and long-range centering interactions (in a synchronous updating framework) 
is not enough for the system to develop directional switching (the alternating flock 
phase); see, e.g., all the cases with = in figures |4| [sj |6| |9] and [Toj Even more, 
the extra stochasticity coming form the probabilistic application of these interactions in 
Version 1 is still not enough. Nonetheless, both versions of our model have built-in the 
capacity to develop directional switching, that gets unleashed as soon as an extra source 
of noise is added to the system. In our case, the simple additive noise term in ([T]) 
and ([s]) does the trick, which leads us to the following question: Why does the noisier 
Version 1 behave qualitatively in the same way as Version 2 regarding the collective 
behaviour of the system, including the directional switching, and the absence of it when 
the additive noise is not present? 

To answer this question, with an eye on the arguments provided above, we propose 
the study of a free flock. For this, we will relax the periodic boundary conditions in both 
versions of our model, allowing the system to freely move in all of position space. We 
can do so, since the cohesion of our system is guaranteed by the long-range character 
of our centering interaction that contrasts with the local interactions considered in the 
previously mentioned models. This also leads us back to the case of the third row in 
figure [9] where, even for the strong noise case (77 = 7 in the figure), the alternating flock 
phase is absent. In clarifying the origin of directional switching in our model, we hope 
to shed some light on the origin of this kind of state in the models discussed above, but 
before, there are some extra features of our model that must be mentioned: 

• First, none of the properties of our model (in both of its versions) depend on the 
average particle density pa, but only on the number of particles N as will be made 
clear in the next subsection. This is a consequence of the long-range character of 
the centering interaction that we have chosen, and resembles the properties of some 
other models with this kind of interactions (see, e.g., [ini[21j) — in particular, the 
Mikhailov and Zanette model that inspired ours. 

• Second, for values of /i > 0.2, both versions of our model are unable to develop an 
alternating flock phase in the steady state, at least for systems with = 1000 as 
shown in figure [llta), where r] = 7 and A = 1. 




This leaves us with the fact that, even with all of the sources of stochasticity active in 
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Figure 11. In (a), (b) and (c) the dynamics of Version 1 and Version 2 of our model, for 
N — 1000, vo ~ 0.1 and /i = 0.2, is presented for different values of the noise intensity 
r] and the alignment interaction radius A. The darker grey level represents higher 
particle density. Periodic boundary conditions have been relaxed, as the particles are 
allowed to freely move in all of position space. They were only applied to the snapshots 
we show, in order to limit them to a finite region of space: in this case, L — 1000. 
The systems were allowed to evolve 30000 time steps. In (a) ?7 = 7 and A = 1. In (b) 
rj — 9 and A = 1. In (c) = 7 and A — 0.8. (d) Semi-log plots of the normalized 
distribution P' {y) of the individual particles distance to the centre of mass, measured 
in the steady state (starting after 10000 time steps of evolution) for the different cases 
(see the text for more details). 



Version 1 and Version 2 or our model as explained above, a strong-enough centering 
interaction can suppress the alternating flock phase. In this way, given a fixed number 
of particles N , we will consider two quantities to vary in order to recover it: the noise 
intensity rj and the interaction radius A. Regarding the latter, let us not forget that the 
condition vq ^ 1 for the small velocity regime must be fulfilled (At = 1 and vq = 0.1 
are kept as before). 

In figure [TT|^b), we have increased the amplitude of the noise from rj = 7 to 9 
compared to the case of figure [lT|^ a). As shown in the figure, the alternating fiock phase 
is recovered in both versions of our model with the same qualitative results. This is also 



true for the case of figure 11 ^c), where we have decreased the interaction radius from 



A = 1 of figure 11 a) to 0.8 — this value is still large enough to keep the system within 
the small velocity regime. At this point, it is illustrative to see what happens with 
the steady-state-accumulated normalized distribution P'{y) of the individual particles 
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distance to the centre of mass, defined as: 

yi{t) = - sign[McMW] ■ [xiit) - xcuit)]- (9) 
The results for the different cases are shown in figure [lT|(d) for both versions of our 



model. One can appreciate that the spreading of particles around the centre of mass 
increases with the noise intensity r], when comparing the black solid and dotted curves 
versus the dashed one. Our results show that the width of the distribution P'{y) also 
increases with the number of particles N, while it decreases as the centering parameter 
/i becomes larger. In consequence, for the cases with A = 1 (cases (a) and (b) in figure 



11 ), as noise increases from 77 = 7 to 9, the alternating flock phase is recovered. On the 



other hand, for the cases with rj = 7 (cases (a) and (c) in figure 11), the width of P'{y) 
for Version 1 and Version 2 are practically the same (black solid and dotted curves), 
regardless of the value of A. Thus, reducing the value of the latter from 1 to 0.8 allows 
us to recover the alternating flock phase. 

From this analysis, one can conclude that the important quantities to compare, 
in both versions of our model, is the interaction radius of the alignment interaction 
A with the width of the distribution P'{y)- Moreover, one can also understand why 
these two versions behave qualitatively in the same way regarding the collective motion 
of the group, even though Version 1 has an extra source of noise. Looking at the 
width of the distribution P'{y) for = and A = 1, filled solid grey curves in figure 
[IT]^d), one can notice that the width of P'{y) is larger for Version 1 than for Version 2. 
However, it is much less than the spread induced by the simple additive noise in ([T]) 
and (|8]) for ?7 > 0. In contrast, the extra source of noise in Version 1, coming from the 
probabilistic application of alignment and centering interactions, has stronger effects on 
the fluctuations of the order parameter 0(t) (see figure |6](d)), i.e., on the values of the 
individual particles velocities Ui as shown by the widths of the peaks in the distribution 
P{u) for the cases presented in figure |6|^c) and discussed before. 

Let us not forget that in all of the models derived from the Czirok et al model 
discussed so far (including ours) that present directional switching, the sources of noise 
are all introduced in the velocity terms of the particles. Moreover, in these models 
(excluding ours), the ranges of the interactions are all local. What has not been made 
clear in previous studies is the impact of the different sources of stochasticity on the 
diffusion of the particles in position space in comparison to the ranges of the different 
interactions considered. This may be important in relation to how correlations in the 
motion of particles are kept in time for different systems and the directional switching 
phenomena. Further studies are being made in this direction. 

Additionally, another difference between both versions of our model becomes 
evident when looking at the order parameter 0ss in figures [2|b) and |2](c), and their 
corresponding details in figures ge) and|2{f), as the curves corresponding to Version 2 
of our model (curves with clear symbols) are always below to the corresponding ones 
for Version 1 (curves with solid black symbols), even for different system sizes. This is 
somewhat counterintuitive in view of the previous results, as the noisier version of our 
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model, Version 1, develops higher mean speeds than Version 2 in general. As explained 
before, the order parameter is always smaller than 1 for both versions of our model, 
even in the absence of the additive noise [r] = 0). A plausible explanation leads us 
back to the two peaks shown by P{u) in the steady state for this case (see figure [6]^c)). 
While those for Version 1 of our model are symmetric, they are also wider and clearly 
asymmetric in height, i.e., there is a peak that dominates the distribution centered 
around a larger value than the narrower and more-balanced-in-height peaks of P{u) for 
Version 2. Considering that these peaks correspond to the distribution of the velocities of 
the individual particles, this fact could account for the difference in the order parameter 
0SS for the different cases of the two versions of our model — in principle, one could 
obtain the accumulated steady state order parameter as (pss = \ uP{u) du\. 

At this point, it is worth to draw our attention to the cases with rj = 1 in the upper 
plots of figures |4](d), |5](d) and [lO|^d), where P{u) becomes almost fiat and symmetric, 
very similar to the case without centering (/i = 0) shown in the upper plot of figure 
[TJ^d). In the latter, this distribution corresponds to the strong-noise-induced disordered 
phase, characterized by a random velocity field, where the velocity of every individual 
particle should average to zero throughout the whole evolution of the system. This is 
also true in the presence of centering (/i > 0) if we were to look just at the distribution 
P{u). However, there is a big difference when we focus on the collective behaviour of 
the system. Whereas the mean velocity of the group tends to zero for /i = 0, as /i 
increases, the mean velocity of the group also increases with the noise amplitude as 
explained before. One could wonder about the diffusion of the centre of mass for /i > 
and the net displacement it may develop in the long run, or at least, the mean square 
displacement. 



3.3. Scaling 

Let us now briefly analyze how the properties of both versions of our model scale with 
the size of the system A^. First, we focus on the accumulated order parameter 0ss- As 
shown in figures [2]^b) and|2]^c), and their corresponding details [2]^e) and[2]^f), for weak 
noise (including = 0) the curves are rather fiat regarding their dependance on rj and 
show almost no difference as A^ changes. This is followed by a crossover region and, as 
7] increases, the order parameter increases with A^. In second place, we have already 



mentioned, in subsection |3.2[ that the spread of the position of particles around the 
centre of mass (the width of the distribution P'{y)) increases with A^ as well. This leads 
to the conclusion that larger systems will tolerate larger values of /i before the alternating 
flock phase is killed, given the interaction radius A of the alignment interaction is kept 
constant. In third place comes the frequency of the changes in direction in the same 
phase. As shown in figure |2]^d) for both versions of our model (solid black lines for 
Version 1 and solid grey lines for Version 2), as A^ increases, the frequency of the changes 
in direction becomes smaller. In other words, the intervals between changes in direction 
increase with the size of the system. This is consistent with the trends reported for 
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some other models that also show directional switching [6l HE] . 

Overall, both versions of our model scale in the same way. Moreover, aside from 
the effects of the extra source of noise in Version 1, and fluctuations intrinsic to 
any numerical analysis, the functional dependence of the quantities mentioned in the 
previous paragraph on some of the parameters (e.g., the dependence of 0ss on 77) seem 
to be quantitatively equivalent across different system sizes. On the other hand, since 
we only analyzed three system sizes, = 500,1000,4000, it is no possible for us to 
provide how these quantities depend explicitly on A^. We will leave this analysis for a 
future work, given that all our results point to the fact that it is equivalent to study 
Version 1 or Version 2 of our model regarding the collective behaviour of the system 
and its scaling. 

4. Conclusions 

We have successfully introduced a long-range centering interaction in the model of Czirok 
et al [I9j via two approaches for applying centering and alignment, that derives in two 
characteristic condensed states. The first, for weak noise, consists in a coherent state 
with broken symmetry. Here, the fluctuations induced by the combination of centering 
and alignment in the mean velocity of the densest cluster seem to be stronger than the 
noise-induced ones. Thus, the system reaches a steady-state speed smaller than the 
prescribed one even in the absence of noise. The second, for strong noise, consists in 
an alternating flock phase with non-vanishing mean velocity. In this case, the centering 
interaction combined with the alignment one provides an ordering mechanism for the 
densest cluster that effectively performs oscillations around the centre of mass as it 
changes its direction of motion alternatively. This phase resembles alternating states 
shown by some other models [5l El lEl HZl [18], but also the oscillatory noisy state of the 
Mikhailov and Zanette model [T9] . 

The two versions of our model, one where alignment and centering interactions are 
applied in a probabilistic way and one deterministic, have allowed us to address the 
directional switching phenomenon that seems to be an intrinsic property of the motion 
of animal groups as recent experimental results show O |6]. Even though the origin 
of the stochasticity in real systems is far from clear, we were able to introduce and 
understand a new mechanism for directional switching in ID SPP systems with long- 
range centering. From our results, new insight has been gained into the effects of different 
sources of stochasticity on the diffusion of particles in position and velocity spaces, and 
its relation with the ranges of the different interactions considered for the development 
of alternating flock phases. In these terms, we have shown how the alternating flock 
phase of our model can be suppressed and how to recover it. 

On the other hand, due to the long-range character of the centering interaction 
introduced, the properties of our model do not depend on the mean density of particles, 
but only on the number of particles. In consequence, periodic boundary conditions, 
typically used in ID SPP models ^ El [HI [13 [SI], can be dropped as the system is 
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able to maintain its cohesion. Nonetheless, the scaling properties of our model show 
the same trends as those reported for some other models [6l [161 HI]- A more detailed 
scaling analysis is in the works in order to determine their functional dependence. 

Finally, as our model does not show a truly disordered phase (with a random 
velocity field), it would be worth to study the diffusion of the centre of mass, e.g., its 
net and mean square displacements throughout the time evolution of the system. We 
believe our results may be relevant for the understanding of the directional switching 
phenomenon and, in general, for the theory of flocking. 
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